This is an R Markdown
Notebook. When you execute code within the notebook, the results appear
beneath the code.
Try executing this chunk by clicking the Run button within
the chunk or by placing your cursor inside it and pressing
Ctrl+Shift+Enter.
###Simulating Trajectories
tmax <- 10
X1_0 <- 600
X2_0 <- 30
X3_0 <-10^5
rho <- 0.108
delta <- 0.5
eta <- 9.5*10^(-6)
lambda <- 36
N1 <- 1000
C <- 3
sigma1 <- 0.1
sigma2 <- 0.1
sigma2 <- 0.1
make_path <- function(tmax, x1_0, x2_0, x3_0, rho, delta, eta, lambda, N1, C, sigma1, sigma2, sigma3, deltat, n){
#Initializing
nosteps <- tmax/deltat
x1 <- numeric(nosteps + 1)
x2 <- numeric(nosteps + 1)
x3 <- numeric(nosteps + 1)
x <- matrix(data = c(x1, x2, x3), ncol = 3)
x[1,] <- c(x1_0, x2_0, x3_0)
sigma <- c(sigma1, sigma2, sigma3)
#Simulationg the driving BM
dW1 <- rnorm(n = nosteps, mean = 0, sd = sqrt(deltat))
dW2 <- rnorm(n = nosteps, mean = 0, sd = sqrt(deltat))
dW3 <- rnorm(n = nosteps, mean = 0, sd = sqrt(deltat))
dW <- matrix(data = c(dW1, dW2, dW3), ncol = 3)
#Defining drift function
drift <- function(x){
x1 <- x[1]
x2 <- x[2]
x3 <- x[3]
y1 <- lambda - rho*x1 - eta*x1*x3
y2 <- eta*x1*x3 - delta*x2
y3 <- N1*delta*x2 - C*x3
y <- c(y1, y2, y3)
y
}
#Iterating
for (k in (1:(nosteps))){
x[k+1,] <- x[k,] + drift(x[k,])*deltat + sigma*x[k,]*dW[k,] + 0.5*sigma*x[k,]*sigma*(dW[k,]*dW[k,] - deltat)
}
#Creating output and adding time variable
t <- (0:nosteps)*deltat
out <- data.frame(x, t)
colnames(out) <- c("x1", "x2", "x3", 't')
out
}
subsample <- function(x){
x <- as.matrix(x)
x_fine <- matrix(data = 0, nrow = 500, ncol = 4)
for (i in (1:500)){
x_fine[i,] <- x[20*(i-1) + 1,]
}
x_fine <- as.data.frame(x_fine)
colnames(x_fine) <- c('x1','x2','x3','t')
x_coarse <- matrix(data = 0, nrow = 100, ncol = 4)
for (i in (1:100)){
x_coarse[i,] <- x[100*(i-1) + 1,]
}
x_coarse <- as.data.frame(x_coarse)
colnames(x_coarse) <- c('x1','x2','x3','t')
out <- list(x_fine, x_coarse)
out
}
path <- make_path(tmax = 10, x1_0 = 600, x2_0 = 30, x3_0 = 10^5, rho = 0.108, delta = 0.5, eta = 9.5*10^(-6), lambda = 36, N1 = 1000, C = 3, sigma1 = 0.1, sigma2 = 0.1, sigma3 = 0.1, deltat = 0.001)
sampled_path <- subsample(path)
x_fine <- sampled_path[[1]]
x_coarse <- sampled_path[[2]]
###Plotting the trajectory
path_plot_fine <- plot_ly(x_fine, x = ~x1, y = ~x2, z = ~x3,
type = 'scatter3d', mode = 'lines+markers',
line = list(width = 6, color = ~t, colorscale = 'Viridis'),
marker = list(size = 3.5, color = ~t, colorscale = 'Greens', cmin = -20, cmax = 50))
path_plot_fine
path_plot_coarse <- plot_ly(x_coarse, x = ~x1, y = ~x2, z = ~x3,
type = 'scatter3d', mode = 'lines+markers',
line = list(width = 6, color = ~t, colorscale = 'Viridis'),
marker = list(size = 3.5, color = ~t, colorscale = 'Greens', cmin = -20, cmax = 50))
path_plot_coarse
Euler-Maruyama Estimator
For making inference, we need to have a likelihood for a ny given
parameter \(\theta\). This will be the
sum of the likelihoods of each step, under the EM-hypothesis. Thus, for
each increment, we need to assign a likelihood of that increment, given
the parameter vector. In this case, the stochastic increments to each
coordinate are independent, due to the diagonal structure of \(\Sigma\). So the stochastic quantities here
are Normal and Chi-square distributed. But we have a sum here, but one
is a deterministic transformation of the other. But there might be more
that one dW that give the same increment. In that case i suppose we have
to assign both the likelihoods. Do we need to take a convolution to find
the density of the sum? Apparently, we are to use a Gaussian
approximation as given on the slides. The books probably cover this as
well. Upon looking at the slides again, the EM-based estimator seems
quite doable.
drift <- function(x){
x1 <- x[1]
x2 <- x[2]
x3 <- x[3]
y1 <- lambda - rho*x1 - eta*x1*x3
y2 <- eta*x1*x3 - delta*x2
y3 <- N1*delta*x2 - C*x3
y <- c(y1, y2, y3)
y
}
likelihood_step <- function(theta, x, x_next){
# need to assign subparameters of theta to the names used in the code
y <- drift(x)
a <- log(((2*pi*deltat)^3) * (x[1]*x[2]*x[3]*sigma1*sigma2*sigma3)^2)
b <- ((x_next - x + y*deltat)^2)*c((x[1]*sigma1)^(-2), (x[2]*sigma2)^(-2), (x[3]*sigma3)^(-2))/deltat
}
a+b
Add a new chunk by clicking the Insert Chunk button on the
toolbar or by pressing Ctrl+Alt+I.
When you save the notebook, an HTML file containing the code and
output will be saved alongside it (click the Preview button or
press Ctrl+Shift+K to preview the HTML file).
The preview shows you a rendered HTML copy of the contents of the
editor. Consequently, unlike Knit, Preview does not
run any R code chunks. Instead, the output of the chunk when it was last
run in the editor is displayed.
LS0tDQp0aXRsZTogIlIgTm90ZWJvb2siDQpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sNCi0tLQ0KDQpUaGlzIGlzIGFuIFtSIE1hcmtkb3duXShodHRwOi8vcm1hcmtkb3duLnJzdHVkaW8uY29tKSBOb3RlYm9vay4gV2hlbiB5b3UgZXhlY3V0ZSBjb2RlIHdpdGhpbiB0aGUgbm90ZWJvb2ssIHRoZSByZXN1bHRzIGFwcGVhciBiZW5lYXRoIHRoZSBjb2RlLiANCg0KVHJ5IGV4ZWN1dGluZyB0aGlzIGNodW5rIGJ5IGNsaWNraW5nIHRoZSAqUnVuKiBidXR0b24gd2l0aGluIHRoZSBjaHVuayBvciBieSBwbGFjaW5nIHlvdXIgY3Vyc29yIGluc2lkZSBpdCBhbmQgcHJlc3NpbmcgKkN0cmwrU2hpZnQrRW50ZXIqLg0KDQojIyNTaW11bGF0aW5nIFRyYWplY3Rvcmllcw0KDQpgYGB7cn0NCnRtYXggPC0gMTANClgxXzAgPC0gNjAwDQpYMl8wIDwtIDMwDQpYM18wIDwtMTBeNQ0KcmhvIDwtIDAuMTA4DQpkZWx0YSA8LSAwLjUNCmV0YSA8LSA5LjUqMTBeKC02KQ0KbGFtYmRhIDwtIDM2DQpOMSA8LSAxMDAwDQpDIDwtIDMNCnNpZ21hMSA8LSAwLjENCnNpZ21hMiA8LSAwLjENCnNpZ21hMiA8LSAwLjENCmBgYA0KDQpgYGB7cn0NCm1ha2VfcGF0aCA8LSBmdW5jdGlvbih0bWF4LCB4MV8wLCB4Ml8wLCB4M18wLCByaG8sIGRlbHRhLCBldGEsIGxhbWJkYSwgTjEsIEMsIHNpZ21hMSwgc2lnbWEyLCBzaWdtYTMsIGRlbHRhdCwgbil7DQogICNJbml0aWFsaXppbmcNCiAgbm9zdGVwcyA8LSB0bWF4L2RlbHRhdA0KICB4MSA8LSBudW1lcmljKG5vc3RlcHMgKyAxKQ0KICB4MiA8LSBudW1lcmljKG5vc3RlcHMgKyAxKQ0KICB4MyA8LSBudW1lcmljKG5vc3RlcHMgKyAxKQ0KICB4IDwtIG1hdHJpeChkYXRhID0gYyh4MSwgeDIsIHgzKSwgbmNvbCA9IDMpDQogIHhbMSxdIDwtIGMoeDFfMCwgeDJfMCwgeDNfMCkNCiAgc2lnbWEgPC0gYyhzaWdtYTEsIHNpZ21hMiwgc2lnbWEzKQ0KICANCiAgI1NpbXVsYXRpb25nIHRoZSBkcml2aW5nIEJNDQogIGRXMSA8LSBybm9ybShuID0gbm9zdGVwcywgbWVhbiA9IDAsIHNkID0gc3FydChkZWx0YXQpKQ0KICBkVzIgPC0gcm5vcm0obiA9IG5vc3RlcHMsIG1lYW4gPSAwLCBzZCA9IHNxcnQoZGVsdGF0KSkNCiAgZFczIDwtIHJub3JtKG4gPSBub3N0ZXBzLCBtZWFuID0gMCwgc2QgPSBzcXJ0KGRlbHRhdCkpDQogIGRXIDwtIG1hdHJpeChkYXRhID0gYyhkVzEsIGRXMiwgZFczKSwgbmNvbCA9IDMpDQogIA0KICAjRGVmaW5pbmcgZHJpZnQgZnVuY3Rpb24NCiAgZHJpZnQgPC0gZnVuY3Rpb24oeCl7DQogICAgeDEgPC0geFsxXQ0KICAgIHgyIDwtIHhbMl0NCiAgICB4MyA8LSB4WzNdDQogICAgeTEgPC0gbGFtYmRhIC0gcmhvKngxIC0gZXRhKngxKngzDQogICAgeTIgPC0gZXRhKngxKngzIC0gZGVsdGEqeDINCiAgICB5MyA8LSBOMSpkZWx0YSp4MiAtIEMqeDMNCiAgICB5IDwtIGMoeTEsIHkyLCB5MykNCiAgICB5DQogIH0NCiAgDQogICNJdGVyYXRpbmcNCiAgZm9yIChrIGluICgxOihub3N0ZXBzKSkpew0KICAgIHhbaysxLF0gPC0geFtrLF0gKyBkcmlmdCh4W2ssXSkqZGVsdGF0ICsgc2lnbWEqeFtrLF0qZFdbayxdICsgMC41KnNpZ21hKnhbayxdKnNpZ21hKihkV1trLF0qZFdbayxdIC0gZGVsdGF0KQ0KICB9DQogIA0KICAjQ3JlYXRpbmcgb3V0cHV0IGFuZCBhZGRpbmcgdGltZSB2YXJpYWJsZQ0KICB0IDwtICgwOm5vc3RlcHMpKmRlbHRhdA0KICBvdXQgPC0gZGF0YS5mcmFtZSh4LCB0KQ0KICBjb2xuYW1lcyhvdXQpIDwtIGMoIngxIiwgIngyIiwgIngzIiwgJ3QnKQ0KICBvdXQNCn0NCmBgYA0KDQpgYGB7cn0NCnN1YnNhbXBsZSA8LSBmdW5jdGlvbih4KXsNCiAgeCA8LSBhcy5tYXRyaXgoeCkNCiAgDQogIHhfZmluZSA8LSBtYXRyaXgoZGF0YSA9IDAsIG5yb3cgPSA1MDAsIG5jb2wgPSA0KQ0KICBmb3IgKGkgaW4gKDE6NTAwKSl7DQogICAgeF9maW5lW2ksXSA8LSB4WzIwKihpLTEpICsgMSxdDQogIH0NCiAgeF9maW5lIDwtIGFzLmRhdGEuZnJhbWUoeF9maW5lKQ0KICBjb2xuYW1lcyh4X2ZpbmUpIDwtIGMoJ3gxJywneDInLCd4MycsJ3QnKQ0KICANCiAgeF9jb2Fyc2UgPC0gbWF0cml4KGRhdGEgPSAwLCBucm93ID0gMTAwLCBuY29sID0gNCkNCiAgZm9yIChpIGluICgxOjEwMCkpew0KICAgIHhfY29hcnNlW2ksXSA8LSB4WzEwMCooaS0xKSArIDEsXQ0KICB9DQogIHhfY29hcnNlIDwtIGFzLmRhdGEuZnJhbWUoeF9jb2Fyc2UpDQogIGNvbG5hbWVzKHhfY29hcnNlKSA8LSBjKCd4MScsJ3gyJywneDMnLCd0JykNCiAgDQogIG91dCA8LSBsaXN0KHhfZmluZSwgeF9jb2Fyc2UpDQogIG91dA0KfQ0KYGBgDQpgYGB7cn0NCnBhdGggPC0gbWFrZV9wYXRoKHRtYXggPSAxMCwgeDFfMCA9IDYwMCwgeDJfMCA9IDMwLCB4M18wID0gMTBeNSwgcmhvID0gMC4xMDgsIGRlbHRhID0gMC41LCBldGEgPSA5LjUqMTBeKC02KSwgbGFtYmRhID0gMzYsIE4xID0gMTAwMCwgQyA9IDMsIHNpZ21hMSA9IDAuMSwgc2lnbWEyID0gMC4xLCBzaWdtYTMgPSAwLjEsIGRlbHRhdCA9IDAuMDAxKQ0Kc2FtcGxlZF9wYXRoIDwtIHN1YnNhbXBsZShwYXRoKQ0KeF9maW5lIDwtIHNhbXBsZWRfcGF0aFtbMV1dDQp4X2NvYXJzZSA8LSBzYW1wbGVkX3BhdGhbWzJdXQ0KYGBgDQoNCiMjI1Bsb3R0aW5nIHRoZSB0cmFqZWN0b3J5DQpgYGB7ciBpbmNsdWRlPUZBTFNFfQ0KbGlicmFyeShwbG90bHkpDQpgYGANCmBgYHtyfQ0KcGF0aF9wbG90X2ZpbmUgPC0gcGxvdF9seSh4X2ZpbmUsIHggPSB+eDEsIHkgPSB+eDIsIHogPSB+eDMsDQogICAgICAgICAgICAgICAgICAgICAgICAgIHR5cGUgPSAnc2NhdHRlcjNkJywgbW9kZSA9ICdsaW5lcyttYXJrZXJzJywNCiAgICAgICAgICAgICAgICAgICAgICAgICAgbGluZSA9IGxpc3Qod2lkdGggPSA2LCBjb2xvciA9IH50LCBjb2xvcnNjYWxlID0gJ1ZpcmlkaXMnKSwNCiAgICAgICAgICAgICAgICAgICAgICAgICAgbWFya2VyID0gbGlzdChzaXplID0gMy41LCBjb2xvciA9IH50LCBjb2xvcnNjYWxlID0gJ0dyZWVucycsIGNtaW4gPSAtMjAsIGNtYXggPSA1MCkpDQpwYXRoX3Bsb3RfZmluZQ0KYGBgDQpgYGB7cn0NCnBhdGhfcGxvdF9jb2Fyc2UgPC0gcGxvdF9seSh4X2NvYXJzZSwgeCA9IH54MSwgeSA9IH54MiwgeiA9IH54MywNCiAgICAgICAgICAgICAgICAgICAgICAgICAgdHlwZSA9ICdzY2F0dGVyM2QnLCBtb2RlID0gJ2xpbmVzK21hcmtlcnMnLA0KICAgICAgICAgICAgICAgICAgICAgICAgICBsaW5lID0gbGlzdCh3aWR0aCA9IDYsIGNvbG9yID0gfnQsIGNvbG9yc2NhbGUgPSAnVmlyaWRpcycpLA0KICAgICAgICAgICAgICAgICAgICAgICAgICBtYXJrZXIgPSBsaXN0KHNpemUgPSAzLjUsIGNvbG9yID0gfnQsIGNvbG9yc2NhbGUgPSAnR3JlZW5zJywgY21pbiA9IC0yMCwgY21heCA9IDUwKSkNCnBhdGhfcGxvdF9jb2Fyc2UNCmBgYA0KDQojIyMgRXVsZXItTWFydXlhbWEgRXN0aW1hdG9yDQpGb3IgbWFraW5nIGluZmVyZW5jZSwgd2UgbmVlZCB0byBoYXZlIGEgbGlrZWxpaG9vZCBmb3IgYSBueSBnaXZlbiBwYXJhbWV0ZXIgJFx0aGV0YSQuIFRoaXMgd2lsbCBiZSB0aGUgc3VtIG9mIHRoZSBsaWtlbGlob29kcyBvZiBlYWNoIHN0ZXAsIHVuZGVyIHRoZSBFTS1oeXBvdGhlc2lzLiBUaHVzLCBmb3IgZWFjaCBpbmNyZW1lbnQsIHdlIG5lZWQgdG8gYXNzaWduIGEgbGlrZWxpaG9vZCBvZiB0aGF0IGluY3JlbWVudCwgZ2l2ZW4gdGhlIHBhcmFtZXRlciB2ZWN0b3IuIEluIHRoaXMgY2FzZSwgdGhlIHN0b2NoYXN0aWMgaW5jcmVtZW50cyB0byBlYWNoIGNvb3JkaW5hdGUgYXJlIGluZGVwZW5kZW50LCBkdWUgdG8gdGhlIGRpYWdvbmFsIHN0cnVjdHVyZSBvZiAkXFNpZ21hJC4gU28gdGhlIHN0b2NoYXN0aWMgcXVhbnRpdGllcyBoZXJlIGFyZSBOb3JtYWwgYW5kIENoaS1zcXVhcmUgZGlzdHJpYnV0ZWQuIEJ1dCB3ZSBoYXZlIGEgc3VtIGhlcmUsIGJ1dCBvbmUgaXMgYSBkZXRlcm1pbmlzdGljIHRyYW5zZm9ybWF0aW9uIG9mIHRoZSBvdGhlci4gQnV0IHRoZXJlIG1pZ2h0IGJlIG1vcmUgdGhhdCBvbmUgZFcgdGhhdCBnaXZlIHRoZSBzYW1lIGluY3JlbWVudC4gSW4gdGhhdCBjYXNlIGkgc3VwcG9zZSB3ZSBoYXZlIHRvIGFzc2lnbiBib3RoIHRoZSBsaWtlbGlob29kcy4gRG8gd2UgbmVlZCB0byB0YWtlIGEgY29udm9sdXRpb24gdG8gZmluZCB0aGUgZGVuc2l0eSBvZiB0aGUgc3VtPyBBcHBhcmVudGx5LCB3ZSBhcmUgdG8gdXNlIGEgR2F1c3NpYW4gYXBwcm94aW1hdGlvbiBhcyBnaXZlbiBvbiB0aGUgc2xpZGVzLiBUaGUgYm9va3MgcHJvYmFibHkgY292ZXIgdGhpcyBhcyB3ZWxsLg0KVXBvbiBsb29raW5nIGF0IHRoZSBzbGlkZXMgYWdhaW4sIHRoZSBFTS1iYXNlZCBlc3RpbWF0b3Igc2VlbXMgcXVpdGUgZG9hYmxlLg0KYGBge3J9DQpkcmlmdCA8LSBmdW5jdGlvbih4KXsNCiAgICB4MSA8LSB4WzFdDQogICAgeDIgPC0geFsyXQ0KICAgIHgzIDwtIHhbM10NCiAgICB5MSA8LSBsYW1iZGEgLSByaG8qeDEgLSBldGEqeDEqeDMNCiAgICB5MiA8LSBldGEqeDEqeDMgLSBkZWx0YSp4Mg0KICAgIHkzIDwtIE4xKmRlbHRhKngyIC0gQyp4Mw0KICAgIHkgPC0gYyh5MSwgeTIsIHkzKQ0KICAgIHkNCiAgfQ0KbGlrZWxpaG9vZF9zdGVwIDwtIGZ1bmN0aW9uKHRoZXRhLCB4LCB4X25leHQpew0KICAjIG5lZWQgdG8gYXNzaWduIHN1YnBhcmFtZXRlcnMgb2YgdGhldGEgdG8gdGhlIG5hbWVzIHVzZWQgaW4gdGhlIGNvZGUNCiAgeSA8LSBkcmlmdCh4KQ0KICBhIDwtIGxvZygoKDIqcGkqZGVsdGF0KV4zKSAqICh4WzFdKnhbMl0qeFszXSpzaWdtYTEqc2lnbWEyKnNpZ21hMyleMikNCiAgYiA8LSAoKHhfbmV4dCAtIHggKyB5KmRlbHRhdCleMikqYygoeFsxXSpzaWdtYTEpXigtMiksICh4WzJdKnNpZ21hMileKC0yKSwgKHhbM10qc2lnbWEzKV4oLTIpKS9kZWx0YXQNCiAgfQ0KICBhK2INCmBgYA0KDQoNCkFkZCBhIG5ldyBjaHVuayBieSBjbGlja2luZyB0aGUgKkluc2VydCBDaHVuayogYnV0dG9uIG9uIHRoZSB0b29sYmFyIG9yIGJ5IHByZXNzaW5nICpDdHJsK0FsdCtJKi4NCg0KV2hlbiB5b3Ugc2F2ZSB0aGUgbm90ZWJvb2ssIGFuIEhUTUwgZmlsZSBjb250YWluaW5nIHRoZSBjb2RlIGFuZCBvdXRwdXQgd2lsbCBiZSBzYXZlZCBhbG9uZ3NpZGUgaXQgKGNsaWNrIHRoZSAqUHJldmlldyogYnV0dG9uIG9yIHByZXNzICpDdHJsK1NoaWZ0K0sqIHRvIHByZXZpZXcgdGhlIEhUTUwgZmlsZSkuDQoNClRoZSBwcmV2aWV3IHNob3dzIHlvdSBhIHJlbmRlcmVkIEhUTUwgY29weSBvZiB0aGUgY29udGVudHMgb2YgdGhlIGVkaXRvci4gQ29uc2VxdWVudGx5LCB1bmxpa2UgKktuaXQqLCAqUHJldmlldyogZG9lcyBub3QgcnVuIGFueSBSIGNvZGUgY2h1bmtzLiBJbnN0ZWFkLCB0aGUgb3V0cHV0IG9mIHRoZSBjaHVuayB3aGVuIGl0IHdhcyBsYXN0IHJ1biBpbiB0aGUgZWRpdG9yIGlzIGRpc3BsYXllZC4NCg==